#!/bin/bash
#
# Developed by Fred Weinhaus 8/25/2009 .......... revised 1/16/2020
#
# ------------------------------------------------------------------------------
# 
# Licensing:
# 
# Copyright © Fred Weinhaus
# 
# My scripts are available free of charge for non-commercial use, ONLY.
# 
# For use of my scripts in commercial (for-profit) environments or 
# non-free applications, please contact me (Fred Weinhaus) for 
# licensing arrangements. My email address is fmw at alink dot net.
# 
# If you: 1) redistribute, 2) incorporate any of these scripts into other 
# free applications or 3) reprogram them in another scripting language, 
# then you must contact me for permission, especially if the result might 
# be used in a commercial or for-profit environment.
# 
# My scripts are also subject, in a subordinate manner, to the ImageMagick 
# license, which can be found at: http://www.imagemagick.org/script/license.php
# 
# ------------------------------------------------------------------------------
# 
####
#
# USAGE: normcrosscorr [-s] [-p] [-m mode] [-c color] [-t type] 
# smallfile largefile corrfile [matchfile]
# USAGE: normcrosscorr [-h or -help]
#
# OPTIONS:
#
# -s                stretch the correlation surface to full dynamic range so 
#                   that the best match is full white; default is unstretched
# -p                apply a pseudocolor to the correlation surface image; 
#                   default is no pseudocoloring
# -m      mode      mode for matchfile output; choices are: draw, overlay or 
#                   best; draw colored box at best match location, overlay
#                   the small image at match location on a one half transparent
#                   large image or output the best match subsection; 
#                   default=draw
# -c      color     color to use for drawing box on large image where best 
#                   matched subsection was found; default=black
# -t      type      type of approach used to combine color correlation results 
#                   into a single grayscale surface image; choices are: gray,
#                   rec709luma, rec601luma, average and rms; default=gray
#
###
#
# NAME: NORMCROSSCORR 
# 
# PURPOSE: To compute the normalized cross correlation surface to find where 
# a small image best matches within a larger image.
# 
# DESCRIPTION: NORMCROSSCORR computes the normalized cross correlation surface 
# (image) to find where a small (first) image best matches within a larger 
# (second) image. Since the matching may differ for each channel, the output  
# correlation image will be converted to grayscale. Any alpha channel on either 
# input image will be removed automatically before processing. Values in the 
# correlation surface can vary between +1 and -1, with a perfect match 
# being +1. If the correlation surface result is saved to an image format 
# that does not support negative values, the correlation surface will be 
# clamped so that all negative values are zero.
# 
# 
# OPTIONS: 
# 
# -s ... Stretch the normalized cross correlation surface image to full 
# dynamic range. Default is no stretch.
# 
# -p ... Apply a pseudocoloring to the normalized cross correlation surface 
# image where red corresponds to the highest values and purple to the lowest 
# values. Default is no pseudocoloring.
# 
# -m mode ... MODE is the layout mode for the optional matchfile image.
# Choices are draw (or d), overlay (or o) or best (or b). Draw simply draws 
# a colored box outline at the best match subsection in the larger image. 
# Overlay inserts the small image at the match location of a 30% opaque 
# version of the larger image. Best outputs the subsection of the larger image 
# that best matches the smaller image. The default="draw". Ignored if no 
# matchfile specified.
# 
# -c color ... COLOR is the color to use to draw the outline of the best 
# matching subsection in the larger image when mode=draw. Any valid IM color 
# specification may be used. The default=black.
# 
# -t type ... TYPE of approach used to combine color correlation results into 
# a single grayscale surface image. The choices are: gray, rec709luma, 
# rec601luma, average and rms. The default=gray. Note that average and rms 
# require IM 6.8.5.5 or higher.
# 
# REQUIREMENTS: IM version 6.5.4-7 or higher, but compiled with HDRI enabled 
# in any quantum level of Q8, Q16 or Q32. Also requires the FFTW delegate 
# library.
# 
# LIMITATIONS: Artifacts can be produced in the correlation surface image 
# when there are constant color regions in the larger image of equal or larger 
# size than the smaller image. This occurs, because the standard deviation in 
# those regions is zero or near zero, thereby causing a divide by (near) zero 
# situation, which in turn cause a larger than 1 correlation value. Some 
# mitigation of this is performed in the script.
# 
# REFERENCES: 
# F. Weinhaus and G. Latshaw, Edge Extraction Based Image Correlation, Proceedings SPIE, Vol. 205, 67 - 75, 1979.
# http://en.wikipedia.org/wiki/Normalized_cross-correlation#Normalized_cross-correlation
# 
# See http://www.fmwconcepts.com/imagemagick/fourier_transforms/fourier.html 
# for more details about the Fourier Transform with ImageMagick.
# 
# CAVEAT: No guarantee that this script will work on all platforms, 
# nor that trapping of inconsistent parameters is complete and 
# foolproof. Use At Your Own Risk. 
# 
# 
######
#

# set default values
stretch="no"		#yes or no
pseudocolor="no"    #yes or no
mode="draw"			#draw or overlay
color="black"		#any valid IM color
type="gray"
transp=0.3			#transparency of large image in mode=overlay

# set directory for temporary files
tmpdir="/tmp"

# set up functions to report Usage and Usage with Description
PROGNAME=`type $0 | awk '{print $3}'`  # search for executable on path
PROGDIR=`dirname $PROGNAME`            # extract directory of program
PROGNAME=`basename $PROGNAME`          # base name of program
usage1() 
	{
	echo >&2 ""
	echo >&2 "$PROGNAME:" "$@"
	sed >&2 -e '1,/^####/d;  /^###/g;  /^#/!q;  s/^#//;  s/^ //;  4,$p' "$PROGDIR/$PROGNAME"
	}
usage2() 
	{
	echo >&2 ""
	echo >&2 "$PROGNAME:" "$@"
	sed >&2 -e '1,/^####/d;  /^######/g;  /^#/!q;  s/^#*//;  s/^ //;  4,$p' "$PROGDIR/$PROGNAME"
	}


# function to report error messages
errMsg()
	{
	echo ""
	echo $1
	echo ""
	usage1
	exit 1
	}

# function crossCorr to compute IFT of complex product of (A*)x(B), 
# where A* is complex conjugate
# A*=a1-ia2; B=b1+ib2
# (A*)x(B)=(a1xb1+a2*b2) + i(a1xb2-a2xb1)
crossCorr()
	{
	img1=$1
	img2=$2
	# note both images contain 2 frames
	if [ "$im_version" -ge "06080701" ]; then
	convert \( $img1 -complex conjugate \) $img2 -complex multiply \
		$fouriernorm +ift "$dir/tmp0.mpc"
	else
		convert $img1 $img2 \
			\( -clone 0 -clone 2 -define compose:clamp=false -compose multiply -composite \) \
			\( -clone 1 -clone 3 -define compose:clamp=false -compose multiply -composite \) \
			\( -clone 4 -clone 5 -define compose:clamp=false -compose plus -composite \) \
			\( -clone 0 -clone 3 -define compose:clamp=false -compose multiply -composite \) \
			\( -clone 1 -clone 2 -define compose:clamp=false -compose multiply -composite \) \
			\( -clone 7 -clone 8 +swap -define compose:clamp=false -compose minus -composite \) \
			-delete 0-5,7,8 $fouriernorm +ift "$dir/tmp0.mpc"
	fi
	}

# function to test for minus at start of value of second part of option 1 or 2
checkMinus()
	{
	test=`echo "$1" | grep -c '^-.*$'`   # returns 1 if match; 0 otherwise
    [ $test -eq 1 ] && errMsg "$errorMsg"
	}

# test for correct number of arguments and get values
if [ $# -eq 0 ]
	then
	# help information
   echo ""
   usage2
   exit 0
elif [ $# -gt 12 ]
	then
	errMsg "--- TOO MANY ARGUMENTS WERE PROVIDED ---"
else
	while [ $# -gt 0 ]
		do
			# get parameter values
			case "$1" in
		  -h|-help)    # help information
					   echo ""
					   usage2
					   exit 0
					   ;;
				-s)    # get stretch
					   stretch="yes"
					   ;;
				-p)    # get pseudocolor
					   pseudocolor="yes"
					   ;;
			   	-m)    # mode
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID MODE SPECIFICATION ---"
					   checkMinus "$1"
					   mode=`echo "$1" | tr "[:upper:]" "[:lower:]"`
					   case "$mode" in
					   		draw|d) mode="draw" ;;
					   		overlay|o) mode="overlay" ;;
					   		best|b) mode="best" ;;
					   		*) errMsg "--- MODE=$mode IS NOT A VALID CHOICE ---" ;;
					   esac
					   ;;
			   	-c)    # get color
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID COLOR SPECIFICATION ---"
					   checkMinus "$1"
					   color="$1"
					   ;;
			   	-t)    # type
					   shift  # to get the next parameter
					   # test if parameter starts with minus sign 
					   errorMsg="--- INVALID TYPE SPECIFICATION ---"
					   checkMinus "$1"
					   type=`echo "$1" | tr "[:upper:]" "[:lower:]"`
					   case "$type" in
					   		gray) ;;
					   		rec709luma) ;;
					   		rec601luma) ;;
					   		average) ;;
					   		rms) ;;
					   		*) errMsg "--- TYPE=$type IS NOT A VALID CHOICE ---" ;;
					   esac
					   ;;
				 -)    # STDIN and end of arguments
					   break
					   ;;
				-*)    # any other - argument
					   errMsg "--- UNKNOWN OPTION ---"
					   ;;
		     	 *)    # end of arguments
					   break
					   ;;
			esac
			shift   # next option
	done
	#
	# get infile, filtfile and outfile
	smallfile="$1"
	largefile="$2"
	corrfile="$3"
	matchfile="$4"
fi

# test that infile provided
[ "$smallfile" = "" ] && errMsg "NO SMALL INPUT FILE SPECIFIED"

# test that filtfile provided
[ "$largefile" = "" ] && errMsg "NO LARGE INPUT FILE SPECIFIED"

# test that outfile provided
[ "$corrfile" = "" ] && errMsg "NO CORRELATION FILE SPECIFIED"



# Setup directory for temporary files
# On exit remove ALL -- the whole directory of temporary images
dir="$tmpdir/$PROGNAME.$$"
trap "rm -rf $dir; exit 0" 0
trap "rm -rf $dir; exit 1" 1 2 3 15
#trap "rm -rf $dir; exit 1" ERR
mkdir "$dir" || {
  echo >&2 "$PROGNAME: Unable to create working dir \"$dir\" -- ABORTING"
  exit 10
}

# read the input image and filter image into the temp files and test validity.
convert -quiet "$smallfile" -alpha off +repage "$dir/tmpA1.mpc" ||
	errMsg "--- FILE $smallfile DOES NOT EXIST OR IS NOT AN ORDINARY FILE, NOT READABLE OR HAS ZERO SIZE  ---"

convert -quiet "$largefile" -alpha off +repage "$dir/tmpA2.mpc" ||
	errMsg "--- FILE $largefile DOES NOT EXIST OR IS NOT AN ORDINARY FILE, NOT READABLE OR HAS ZERO SIZE  ---"


# test for valid version of IM
im_version=`convert -list configure | \
	sed '/^LIB_VERSION_NUMBER */!d;  s//,/;  s/,/,0/g;  s/,0*\([0-9][0-9]\)/\1/g' | head -n 1`
[ "$im_version" -lt "06050407" ] && errMsg "--- REQUIRES IM VERSION 6.5.4-7 OR HIGHER ---"
#echo "im_version=$im_version;"

# test for hdri enabled
if [ "$im_version" -ge "07000000" ]; then
	hdri_on=`convert -version | grep "HDRI"`	
else
	hdri_on=`convert -list configure | grep "enable-hdri"`
fi
[ "$hdri_on" = "" ] && errMsg "--- REQUIRES HDRI ENABLED IN IM COMPILE ---"

if [ "$im_version" -ge "07000000" ]; then
	identifying="magick identify"
else
	identifying="identify"
fi

# get image dimensions to be sure that infile1 is smaller than infile2
ws=`$identifying -ping -format "%w" "$dir/tmpA1.mpc"`
hs=`$identifying -ping -format "%h" "$dir/tmpA1.mpc"`
wlo=`$identifying -ping -format "%w" "$dir/tmpA2.mpc"`
hlo=`$identifying -ping -format "%h" "$dir/tmpA2.mpc"`

[ $ws -gt $wlo ] && errMsg "--- SECOND IMAGE MUST BE WIDER THAN FIRST IMAGE ---"
[ $hs -gt $hlo ] && errMsg "--- SECOND IMAGE MUST BE TALLER THAN FIRST IMAGE ---"

# get diff sizes for final crop to avoid best match with wrap around
wd=$((wlo-ws))
hd=$((hlo-hs))

# compute even dimensions for large image
wl=`convert xc: -format "%[fx:2*ceil($wlo/2)]" info:`
hl=`convert xc: -format "%[fx:2*ceil($hlo/2)]" info:`

# test if large image is odd, and if so, then pad to even with the mean
test1=`convert xc: -format "%[fx:(($wl==$wlo) && ($hl==$hlo))?1:0]" info:`
if [ $test1 -eq 0 ]; then
	meanl=`convert "$dir/tmpA2.mpc" -format "%[fx:100*mean.r]\%,%[fx:100*mean.g]\%,%[fx:100*mean.b]\%" info:`
	convert "$dir/tmpA2.mpc" -background "gray($meanl)" -extent ${wl}x${hl} "$dir/tmpA2.mpc"
fi
#echo "test1=$test1; test2=$test2; ws=$ws; hs=$hs; wlo=$wlo; hlo=$hlo; wl=$wl; hl=$hl;"

: '
C = ((S-Ms) X (L-Ml))/(sigmaS*sqrt(Ns(U X L^2) - (U X L)^2)
Note above does not work when try to compensate for std of Large image equal 0, so use following
C = (((S-Ms) X (L-Ml))/Ns) / (sigmaS*sqrt((U X L^2)/Ns - (U X L)^2/(Ns*Ns))
where
A X B = I(F(A*)F(B)] and A* is complex conjugate of A, F=FFT and I=IFT
L is large image.
L-Ml is the mean subtracted large image. 
S-Ms is the mean subtracted small image. 
U is a unit image (value=1) 
Both S-Ms and U are padded at right and bottom to size of L. 

Note numerator can be expanded from spatial format to show that it is
sum[(S-Ms)(L)] - Ml*sum(S-Ms)
But the second term is zero since sum(S-Ms)=sum(S)-N*Ms and both are just the total pixel values in the small image, so
C = (((S-Ms) X (L))/Ns) / (sigmaS*sqrt((U X L^2)/Ns - (U X L)^2/(Ns*Ns))
'

# compute N=wsxhs = total pixels in small image
Ns=`echo "$ws*$hs" | bc`
Nl=`echo "$wl*$hl" | bc`


# get quantumrange
qrange=`convert xc: -precision 10 -format "%[fx:quantumrange]" info:`
qrangesqr=`echo "scale=10; sqrt($qrange)" | bc`


# get factors
fact1=`echo "scale=10; 1/$Ns" | bc`
fact2=`echo "scale=10; $qrange/($Ns*$Ns)" | bc`

#echo "Ns=$Ns; Nl=$Nl fact1=$fact1; fact2=$fact2; qrange=$qrange; qrangesqr=$qrangesqr;"


# get mean and std of small image
mean=`convert "$dir/tmpA1.mpc" -format "%[fx:100*mean.r]\%,%[fx:100*mean.g]\%,%[fx:100*mean.b]\%" info:`
std=`convert "$dir/tmpA1.mpc" -format "%[fx:100*standard_deviation.r]\%,%[fx:100*standard_deviation.g]\%,%[fx:100*standard_deviation.b]\%" info:`


# ncc has range -1 to +1, so need to scale that to range 0 to quantumrange
# thus we scale by dividing std by quantumrange and 
# note that negatives are clipped by PNG output.
if [ "$im_version" -lt "06050410" ]; then
	#HDRI was unscaled by quantumrange
	# so need the extra factor of quantumrange?
	# not sure why a perfect match is considerably less that quantumrange?
	std=`convert xc: -format "%[fx:$std/(quantumrange)]" info:`
fi
#echo "mean=$mean; std=$std;"

if [ "$im_version" -ge "06080610" ]; then
	fouriernorm="-define fourier:normalize=inverse"
	normfact=""
else
	# default is forward normalization
	fouriernorm=""
	normfact="-evaluate multiply $Nl"
fi
#echo "fouriernorm=$fouriernorm; normfact=$normfact"

# get square of large image and take FFT
# correct for IM normalization by multiplying by quantumrange
convert "$dir/tmpA2.mpc" "$dir/tmpA2.mpc" -define compose:clamp=false -compose multiply -composite -evaluate multiply $qrange $fouriernorm +fft "$dir/tmpL2.mpc"

# take FFT of large image
convert "$dir/tmpA2.mpc" $fouriernorm +fft "$dir/tmpL.mpc"

# subtract mean from small image and pad and take FFT
# unnormalize FFT by multiplying by Nl if normalization is forward
convert "$dir/tmpA1.mpc" \( -size ${ws}x${hs} xc:"rgb($mean)" \) +swap -define compose:clamp=false -compose minus -composite \
	-background black -extent ${wl}x${hl} $fouriernorm +fft $normfact "$dir/tmpS.mpc"

# create identity U image (value=1) and pad and take FFT
# unnormalize FFT by multiplying by Nl if normalization is forward
convert -size ${ws}x${hs} xc:white -background black -extent ${wl}x${hl} $fouriernorm +fft $normfact "$dir/tmpU.mpc"

# create ((S-Ms) X (L))/Ns
crossCorr  "$dir/tmpS.mpc"  "$dir/tmpL.mpc"
convert "$dir/tmp0.mpc" -evaluate multiply $fact1 "$dir/tmpSL.mpc"

# create (U X L^2)/Ns
crossCorr "$dir/tmpU.mpc" "$dir/tmpL2.mpc"
convert "$dir/tmp0.mpc" -evaluate multiply $fact1 "$dir/tmpL2.mpc"

#create (U X L)^2/Ns^2
# multiply by qrange to account for IM normalization in multiply
# so final multiply is qrange/(Ns*Ns)
crossCorr  "$dir/tmpU.mpc"  "$dir/tmpL.mpc"
convert "$dir/tmp0.mpc" "$dir/tmp0.mpc" -define compose:clamp=false -compose multiply -composite -evaluate multiply $fact2 "$dir/tmpL.mpc"

# compute the std of L denominator image
# IM normalizes sqrt in range 0 to 1, then multiplies by qrange
# thus need to divide by sqrt(qrange) to compensate
# make white image
# trap values for std of large image if too close to 0 (when sub region is too constant a color), so no divide by zero
# replacing the use of -fx "(abs(u)<0.002)?1:u"
# note do not need abs since std should be non-negative
convert \( "$dir/tmpL2.mpc" "$dir/tmpL.mpc" +swap -define compose:clamp=false -compose minus -composite \
		-evaluate pow 0.5 -evaluate divide $qrangesqr \) \
	\( -size ${wl}x${hl} xc:white \) \
	\( -clone 0 -threshold 2% -negate \) \
	-compose over -composite "$dir/tmpL.mpc"


# set up grayproc
[ "$type" = "average" -a "$im_version" -lt "06080505" ] && errMsg "--- TYPE=AVERAGE REQUIRES IM 6.8.5.5 OR HIGHER ---"
[ "$type" = "rms" -a "$im_version" -lt "06080505" ] && errMsg "--- TYPE=RMS REQUIRES IM 6.8.5.5 OR HIGHER ---"

if [ "$type" = "gray" ]; then
	gproc="-colorspace gray"
elif [ "$type" = "rec709luma" ]; then
	gproc="-colorspace rec709luma"
elif [ "$type" = "rec601luma" ]; then
	gproc="-colorspace rec601luma"
elif [ "$type" = "average" ]; then
	gproc="-grayscale average"
elif [ "$type" = "rms" ]; then
	gproc="-grayscale rms"
fi


#evaluate normalize cross correlation image
convert "$dir/tmpSL.mpc"  \
	\( "$dir/tmpL.mpc" \( -size ${wl}x${hl} xc:"rgb($std)" \) -define compose:clamp=false -compose multiply -composite \) \
	+swap -compose over -define compose:clamp=false -compose divide -composite -crop ${wd}x${hd}+0+0 +repage $gproc \
	"$dir/tmp0.mpc"



# setup pseudocolor lut
if [ "$pseudocolor" = "yes" ]; then
	convert xc:blueviolet xc:blue xc:cyan xc:green1 \
		xc:yellow xc:orange xc:red +append \
		-filter cubic -resize 256x1 "$dir/tmpP.mpc"
	colorize="$dir/tmpP.mpc -clut"
else
colorize=""
fi

# setup stretch
if [ "$stretch" = "yes" ]; then
	convert \( "$dir/tmp0.mpc" -auto-level \) $colorize "$corrfile"
else
	convert "$dir/tmp0.mpc" $colorize "$corrfile"
fi


#echo "get match"
if [ "$im_version" -ge "06080610" ]; then
	max=`convert $dir/tmp0.mpc -format "%[fx:maxima]" info:`
	coords=`$identifying -define identify:locate=maximum -define identify:limit=1 $dir/tmp0.mpc | tail -n 1 | sed -n 's/^.* \([0-9]*,[0-9]*\)/\1/p'`
	echo ""
	echo "Match Coords: ($coords) And Score In Range 0 to 1: ($max)"
	echo ""
else
	max=`convert "$dir/tmp0.mpc" -format "%[fx:maxima]" info:`
	str=`convert "$dir/tmp0.mpc" -fx "u>=($max-quantumscale)?debug(u):0" null: 2>&1`
	coords=`echo "$str" | sed -n 's/^.*\[\([0-9]*,[0-9]*\)\].*$/\1/p'`
	echo ""
	echo "Match Coords: ($coords) And Score In Range 0 to 1: ($max)"
	echo ""
fi

# compute subsection
ulx=`echo $coords | cut -d,  -f 1`
uly=`echo $coords | cut -d,  -f 2`
subsection="${ws}x${hs}+$ulx+$uly"
#echo "$subsection"


if [ "$matchfile" != "" -a "$mode" = "draw" ]; then
	lrx=$(($ulx+$ws))
	lry=$(($uly+$hs))
#echo "ulx=$ulx; uly=$uly; lrx=$lrx; lry=$lry"
	convert $dir/tmpA2.mpc[${wlo}x${hlo}+0+0] -fill none -stroke "$color" \
		-draw "rectangle $ulx,$uly $lrx,$lry" "$matchfile"
elif [ "$matchfile" != "" -a "$mode" = "overlay" ]; then
	convert \( "$dir/tmpA2.mpc" -alpha on -channel a -evaluate set $transp% +channel \) "$dir/tmpA1.mpc" \
		-geometry "$subsection" -compose over -composite "$matchfile"
elif [ "$matchfile" != "" -a "$mode" = "best" ]; then
	convert $dir/tmpA2.mpc[${ws}x${hs}+${ulx}+${uly}] "$matchfile"
fi
exit 0



